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ABSTRACT 

We study the stability of proto-planetary disks with vertical velocity gradi- 
ents in their equilibrium rotation rates; such gradients are expected to develop 
when dust settles into the midplane. Using a linear stability analysis of a simple 
three-layer model, we show that the onset of instability occurs at a larger value 
of the Richardson number, and therefore for a thicker layer, when the effects 
of Coriolis forces are included. This analysis also shows that even-symmetry 
(midplane-crossing) modes develop faster than odd-symmetry ones. These con- 
clusions are corroborated by a large number of nonlinear numerical simulations 
with two different parameterized prescriptions for the initial (continuous) dust 
distributions. Based on these numerical experiments, the Richardson number 
required for marginal stability is more than an order of magnitude larger than 
the traditional 1/4 value. The dominant modes that grow have horizontal wave- 
lengths of several initial dust scale heights, and in nonlinear stages mix solids 
fairly homogeneously over a comparable vertical range. We conclude that gravi- 
tational instability may be more difficult to achieve than previously thought, and 
that the vertical distribution of matter within the dust layer is likely globally, 
rather than locally, determined. 



Subject headings: hydrodynamics — instabilities — turbulence — planets and 
satellites: formation — planetary systems: protoplanetary disks 
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1. INTRODUCTION 

In the generally accepted model for the formation of planets, the circumstellar disk - 
consisting of material remaining from the formation of the central star - progresses through 
a series of dynamical stages. Over the course of this evolution, the solids contained in the 
disk are believed to settle towards the midplane and form planetesimals, which in turn merge 
to assemble terrestrial planets and cores of gaseous giants. Unfortunately, some of the most 
basic elements of this model are not understood, with the detailed physical processes involved 
at many stages only barely explored. 

One such poorly-understood stage is the formation of planetesimals from interstellar 
dust.^ It is generally accepted that micron-sized particles collide and stick together to form 
larger particles, but it is not clear that this stochastic coagulation process can continue up 
the ladder to generate km-sized planetesimals. Among other difficulties with this model is 
the problem that bodies with approximate sizes between 1 cm and 1 m would drift radially 
inward very fast, as a consequence of drag with the partially pressure-supported gas disk in 
which they are embeded (Weidenschilling 1977). As a way of solving this problem, authors 
have looked for a process that allows particles to grow through this range of sizes in time 
scales similar or smaller to the orbital time (Goldreich & Ward 1973; Weidenschilling 1980; 
Nakagawa, Nakazawa & Hayashi 1981; Nakagawa, Sekiya & Hayashi 1986; WeidenschiUing & 
Cuzzi 1993). For example, Cuzzi, Dobrovolskis & Champney (1993) and Champney, Dobro- 
volskis & Cuzzi (1995) noted that because particles of different sizes are affected differently 
by gas drag they would drift at different speeds, thus increasing their collision rate and 
possibly increasing the particle growth rate. Nevertheless, it remains questionable whether 
solid state forces can grow particles through sticking at the collision velocities likely present 
in proto-planetary disks (see discussion in Youdin & Shu 2002). 

Gravitational instability (GI) has been proposed as an alternative (Safranov 1969; Gol- 
dreich & Ward 1973; Coradini, Federico & Magni 1981; Sekiya 1983). In this model, as 
the dust settles vertically to create a dense midplane subdisk, the spatial particle density 
may increase until the layer crosses the GI density threshold. This model has the advantage 
that particle growth from mm to km sizes happens in a short time scale. Still, there are 
uncertainties in this GI model. For example, several authors (WeidenschiUing & Cuzzi 1993; 
Cuzzi, Dobrovolskis & Hogan 1994; Weidenschilling 1995; Tanga et al. 2004) have pointed 
out that ~ 10 m sized bodies, an intermediate step of the two-staged process proposed by 
Goldreich & Ward (1973), reach the disk midplane with random velocities too large for GI 



^In the outer parts of the disk, gaseous water, ammonia, and other molecules condense into ice, augmenting 
the dust with additional solids. 
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to proceed. 

A more frequently discussed problem with the GI hypothesis arises from the fact that, 
as the sohd particle density increases through vertical settling, the gas/dust mixture be- 
comes less sensitive to large scale gas pressure gradients that are likely to be present in the 
disk. With the loss of radially-supporting pressure gradients, the dust-rich mixture near the 
midplane would orbit the central star with an orbital frequency closer to the keplerian rate, 
and thus develop a velocity shear with its dust-poor envelope. This velocity shear might 
trigger the Kelvin- Helmholtz instability (KHI), which could mix the dust and gas layers and 
stop particle settling before GI sets in (Goldrcich & Ward 1973; Wcidcnschilling 1980). The 
question posed at this point is: what are the appropriate circumstances for the development 
of the KHI in circumstellar disks? Or conversely, is it possible to have a KH-stable disk that 
is susceptible to GI, and so able to form planetesimals within a very short timescale? 

The most favorable circumstance for GI is if there is no turbulence - and no mixing 
of solids - in the absense of KHI. Sekiya & Ishitsu (2000, 2001) performed linear stability 
analyses of vertically shearing disks and concluded that GI can be achieved, avoiding the 
KHI induced mixing, only if the dust-to-gas surface density ratio is much larger than the 
solar abundance. Other authors have reached a similar conclusion while considering the 
potentially dc-stabilizing influence of the Coriolis and tidal forces (Ishitsu & Sekiya 2002, 

2003) , stratification of the medium, two-fiuid effects and radiative cooling (Garaud & Lin 

2004) . Youdin & Shu (2002) further argued that the gas involved in the KHI should be 
able to drag and mix only a finite amount of solids, so that there might be a decoupling 
and precipitation of solids once the dust space density is large enough, even for particle 
sizes at which a good velocity coupling is expected. The precipitated particles would be 
free to continue settling and undergo GI. If a strong central enhancement of dust develops 
(Sekiya 1998) and precipitates out an independent solid layer, then potentially this reduces 
the needed dust-to-gas surface density enhancement required for planetesimals to form via 
GI. Even with the modification proposed by Youdin & Shu (2002), however, the dust-to-gas 
surface density ratio must be enhanced above solar-abundance values ~ 0.01 by an order of 
magnitude. 

Of course, this increased surface density ratio (which can also help accelerate stochastic 
coagulation) does not have to be global. Particle drift in the radial direction (Cuzzi, Dobro- 
volskis & Champney 1993; Youdin & Shu 2002; Haghighipour & Boss 2003; Wcidcnschilling 
2003; Youdin & Chiang 2004), and self-gravitating spiral modes (Rice et al. 2004) could 
provide localized surface density concentrations of dust. Provided that these enhancements 
are long-lived, this could lead to localized GI. 

As discussed above, when dust settles toward the midplane, both (de-stabilizing) shear 
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and (stabilizing) buoyancy increase. In the situation where dust and gas are well-enough 
coupled, they may be considered a single fluid. For a stratified flow, competition between the 
opposing shear and buoyancy effects is traditionally described by the Richardson number. 



where v{z) is the equilibrium shear velocity in the plane perpendicular to e^- Miles (1961) 
and Howard (1961) proved (for incompressible flow in the Boussinesque approximation) that 
Ri < 1/4 somewhere in the layer is a necessary (but not sufficient) condition for instability 
(see also Drazin & Reid 1981, p. 325; a review of the proof is presented in the Appendix 
B of Li, Goodman & Narayan 2003). As Garaud & Lin (2004) argue, however, some of the 
assumptions in this proof, while acceptable in other physical regimes, are not applicable to 
the problem at hand. 

In this paper, we explore the inffuence of an effect not considered in arriving at the 
"Ri > 1/4" stability criterion, namely the Coriolis force. In §2 we present the basic equations 
of the model problem we have defined. In §3 we perform a linear stability analysis of a 
discrete three-layer configuration, consisting of a dust-rich layer surrounded by dust-free 
gas. This analysis relaxes some of the assumptions made by previous work, compares modes 
of instability with even and odd symmetry, and assesses the effect of the Coriohs force. In §4 
we present the results of our numerical experiments, focusing on the effect the Coriolis force 
has on the stability of model disks with two different (continuous) initial dust distributions. 
Finally, in §5 we present our conclusions. 



Consider the reference frame of a fiuid orbiting a central star at a distance Rq, with 
angular velocity ftp- For adiabatic fiow of a gamma-law gas, the equations of hydrodynamics 
in this frame read 



Ri = 



g d\np/dz 



(1) 




2. BASIC EQUATIONS 
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All the symbols are standard. We orient the coordinate system with parallel to ftp, 
pointing radially, and e.y pointing along the azimuthal background flow. Wc shall treat a 
small domain compared to i?o, and disregard the curvature terms in equations (2)- (5). 

In addition to gas with density p^, the flow contains a component of solids (dust + ice). 
We asssume that the gas and solids are strongly coupled, %0Vd — Vg — v. In the strongly 
coupled limit, the dust component of the mixture adds inertia to the flow so that p = pg + pd 
in equations (2) and (3) but the solids do not contribute to the pressure (P = c^p^, where 
Cs is the isothermal sound speed) and energy density in equations (4) and (5). The dust 
component also separately obeys a continuity equation 

^ + V • {pav) = 0. (6) 
For an axisymmetric equilibrium, the x and z components of the momentum equation 

read 



Rn 



Pd{z)+Pg{z] 

1 dPp 

pd{z)+pg{z) dz 



-n\z 




0. 



(7) 
(8) 



Here, d^/dR where VIk is the keplerian orbital frequency, and we approximate 

d^/dz ^ nj^z. 

Prom equation (7), the orbital velocity of the flow in equilibrium depends on the dust 
abundance. In order to further specify our coordinate system, let us deflne as the inertial- 
frame orbital frequency of dust- free gas with a reference value of density at 2; = of pgo- 
With this deflnition of ftp, voy — when Pd — and we can relate to the physical 
parameters of the system. 



1 + 



dPoldR\ 



nl/2 



^K^-oPgO 



(9) 



Here, Pq is the pressure in the dust-free case. 

Solving equation (7) for voy and taking \voy\ <^ ^fRq, we obtain 



P 



1+P 



(10) 
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where ii{z) — pd{z)/pg{z) is the dust abundance, and 

_ dPo/dR\RQ^^=Q 

is the maximum speed that the dust-laden gas can attain, i. e. the difference between the 
keplerian velocity and Q,fRo-^ 

In order to define Pq{R) we will adopt the parametrization of Youdin & Shu (2002), 



E,(i?) = 1700gcm-V.(Ylu) ' ^^^^ 
T{R) = 280°K^(^^) ' (13) 

where is the gas surface density, and T is the temperature of the gas. For the minimum 
solar nebula (MSN) model of Hayashi (1981), fg = fr = ^, P = 3/2 and q = 1/2. In this 
case, VQ^raax/ Cg depends only weakly on i?, 

= 7.3 X 10-2 ( -^T-^i ■ (14) 

Consider now the vertical hydrostatic equilibrium. For vertically-isothermal conditions, 
equation (8) can be written as 

^ = -5|i|l+,Wl. (15) 



S 



which can be solved (analytically or numerically) for Pg{z), given p{z). The total density is 
then p — pg{l+ii). In §4 we shall define initial conditions either by setting ii{z) to a specified 
function, or by requiring the solution p{z) to satisfy certain additional assumptions. 

In describing the equilibrium, we have introduced a natural set of units, namely 
for time, for velocity, and pgo for density. For future reference, we note that in the dust-free 
case, the solution of equation (15) is a gaussian with gas column density E^ = y/2npgoHg, 
where Hg — Cs/VLk- In the MSN model, the gas disk aspect ratio is 



2 



vo,max relates to the rjK parameter used in previous work by fo,max 



3. STABILITY ANALYSIS FOR A THREE-LAYER MODEL 



3.1. Perturbation Equations and General Solution 

Consider a perturbation to tlie above equilibrium witli p = pq + pi,v = Vq + Vi, 
and P — Pq + -Pi- We assume that the flow is incompressible (V • Vi = 0), hnearize 
equations (2)- (5), and, for simplicity in this flrst study, we drop terms in d/dx. Then, 
taking pi, Pi,vix,viy,viz oc exp[i(A;j^y — cut)], the resulting equations can be combined to 
obtain a governing equation for vi^: 



+ 



k, 



LU - kyVQy 



+ <{ dzlnpo 
dzlnpo dzVoy + 



kii 



{UJ - kyVOyY 

dzVOy 



dlviz 



u - kyVoy 



d^voy 



(U - kyVOyY 

ky9iz)dz\npo 



dzVlz 



Viz 



0, (17) 



where dz = d/dz. 



The usual step at this point is to invoke the Boussinesque approximation, dropping 
dzPo everywhere but in the bouyancy term. Because shear and density gradients are directly 
associated for the problem at hand, however, we shall not make this approximation in this 
work (see also discussion in Garaud & Lin 2004). 

Consider now the three layer distribution shown in Figure 1. The middle layer is com- 
posed of the dust and gas mixture, while the top and bottom layers are dust-free. We take 
the initial dust density pd as constant in the middle layer, and the gas density Pg constant 
across all layers (provided is small compared to the total disk thickness Hg, the latter 
is true more generally). This three- layer model is a fair approximation to a more realistic 
continuous distribution for wavelengths larger than the vertical velocity and density gradient 
scales. 

In the chosen reference frame, the top and bottom layers are at rest, while the central 
layer has a constant velocity V = VQ^maxPd/ (Pg+Pd) = "fo.max/^/ (l+A*)- Since d^po = dzVoy = 
within each layer, equation (17) simplifles to 
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Fig. 1. — Schematic of the discrete three layer model. A central gas-dust layer, 
with density pg + pd and thickness 2Hci-, is sandwiched between two semi-infinite 
gaseous atmosphere layers with densities pg. In the reference frame rotating with 
the atmosphere at $7^, the central layer has velocity V at radius Rq. We consider a 
small azimuthal-vertical section of the disk, as shown. 



The solution is 



where 



^ Viz \z\ > Hd, 



dzViz = { ^' ^2 (18) 




Viz = { ^2+6^2^ + yl2_e-^2^ \z\ < Hd, (19) 

z < -Hd, 



^ ky 

VI - 4^^0.2 ^ ^ 

= , (21) 

and where we require that the sign of the square root in equation (20) is to be chosen such 
that its real part ^{Ki) > 0. In order to obtain equation (19), we apply the boundary 
condition vi^ — > when \z\ — > oo. In general, solutions that are either midplane-symmetric 
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or midplane- antisymmetric are possible. Some previous work (Ishitsu & Sekiya 2002, for 
example), has explored only the case in which viz is an odd function of z, so Viz{0) = 0. We 
shall focus much of our attention on the case in which viz is an even function, and explore 
the odd case only for reference. 



3.2. Solutions With vu Even in z 

With viz{z) — viz{—z), A\ — A-i and ^2+ = ^2- — A^. To relate A\ and ^2, we use 
the fact that the vertical displacement of a lagrangian point, 8z = iviz/{uj — kyV), must be 
a continuous function of z across any interface in the flow. Applying this condition at the 
layer discontinuity, this implies 

e-K^H,^^ ^ ^ AJ^K,H, ^ ^-K,H,y (22) 

oj — kyV 

By integrating equation (17) across the z — interface, and substituting the solutions 
(19), we obtain the dispersion relation 



{l + f,){u;-kyVy 



Here, g = g{Hd) = Q%Hd. 



(a; - kyVy 



1/2 



/ 4Q2 \ 1/2 
tanh(X2//d) + 1 - ^ j = kyQi^. (23) 



If we temporarily disregard the Coriolis terms (setting = 0; wc shall discuss the case 
with fii? 7^ in §3.5), then Ki = K2 = ky, and wc can analytically solve equation (23) for 
cu. The condition for instability [that the imaginary part of cu, '^{ur), is non-zero] is 



1 



+ 1 + 1^ 



ta,nh.(kyH(],) 

(see also Drazin & Reid 1981, p. 29). This can be written in alternate form as 



(24) 



kyHd > 



1 + H f 



ta.nh.(kyHd) 



(25) 



We now define an effective Richardson number for the discontinuous distribution in the 
three-layer model. 
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p {v/i^zY 

1 + /X / HdX^ ( Cs 



II \HgJ \VQ 



\,max 



TT 1 + /J / \ / C, 



2 fjL^ V^s/ 
where = 2pdH(i is the dust column density, so that 



TT Hg Srf 



max 



(26) 



With this definition, the criterion for instabihty may be written 

1 + (1 + /.) i^Yih{kyHa) ^ ' 

Notice that, loi p':^ 1 and a moderate kyH^ value, equation (28) reduces to Rigj^ < {kyH^) / p, 
i. e. the layer will still remain stable to wavelengths comparable to Ha even when Rigj^ is 
quite small, for sufficient dust concentration. 



3.3. Solutions With Odd in z 

In the case with Viz{z) — —Viz{z), Ai — —A^ and A2+ — —A2- — A2, and equation (22) 
becomes 



e-Kii/,^^ = ^ ^^(e-^^^'' - e-^^^''). (29) 

UJ rvy V 

The resulting dispersion relation for the odd-symmetry mode is 



{l-rp){u-kyVf 



1 - 



40| 



- kyV)\ 



1/2 



coth(X2i^d) + a;M 1 - 



a;" 



2 \ 1/2 



Again, if we consider the Vlp — ^ case we can obtain criteria for instability. These are the 
same as equations (24), (25), and (28) with coi\i{kyHd) substituted for ta,n\i{kyHd). Notice 
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that both the even and odd modes have the same instabihty condition for short wavelengths 
Rieff{2 + ii)<kyHd or + /^) ^ ^^^2^ ^3^^ 

J- \ 

which is the same as the standard criterion (Chandrasekhar 1981, eq. XI.34, for example) 
for KHI at an interface between two semi-infinite layers. 



3.4. Comparing Even and Odd Modes With Vtp ^Q. 

In order to compare with more realistic models that have a continuous distribution of 
dust, we focus on wavelengths that are comparable to the layer thickness; shorter wavelengths 
will not lead to large scale mixing of the dust layer, and longer wavelengths have smaller 
growth rates. Figure 2 shows the minimum stable layer thickness for three different values 
of kyHd, for a range of dust-to-gas surface density ratios and VQ^rnax/cg = 0.1. Both even- 
and odd-symmetry solutions are shown. In a real system, the dust will (at least initially) 
likely settle toward the midplane at a higher rate than the surface density ratio changes. 
So, we can imagine the disk evolving down along a vertical line in Figure 2, until it reaches 
the stability edge at a relevant wavelength. On a longer time scale, the disk may then move 
along that edge towards larger S^/Sg values as gas is lost to photoevaporation or the dust 
component is increased due to radial drift, for example. 

The slopes of the critical curves at fixed kyH,! can be understood physically. At 
small T,d/Tig, an increase in 'Ld/^g increases /x, thus increasing the velocity shear since 
y — VQ^maxfJ^/ (1 + A*)- This forces an increase in Ha/ Hg in order to keep the layer marginally 
stable; as a consequence, from equation (25), the critical (Hd/Hg) oc (Ed/Eg)-*^/^ at small 
II. For large E^^/Eg, /i ^ 1 and the velocity difference saturates to Vo^max] in this case, the 
dust layer must get thinner to remain marginally stable with increasing dust abundance. At 
large /i, the critical curve therefore follows {Hd/Hg) oc (E^/Eg)"^. These small- and large-/x 
scalings of {Hd/ Hg)crit with E^/Ec, are the same as have been identified by Garaud & Lin 
(2004) in their equations (C3) and (C9), respectively. 

Lines of constant Rie^f are also plotted in Figure 2. Recall that for a continuous dust 
distribution, Ri < 1/4 is a necessary (but not sufficient) condition for the development of 
the KH instability. For our discrete three-layer distribution, Rieif = const, is a good proxy 
for the stability limit of both even and odd modes at fixed kyHd, but only at low E^/E^. 
For kyHd = n/4:, we find that the critical Rigff ~ 0.3 at ^d/^g ^ 0.01, close to the Ri = 1/4 
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Fig. 2. — Minimum stable dust layer thickness for global modes with ftp = (i. e. 
neglecting Coriolis forces), as a function of the dust-to-gas surface density ratio. 
Both even (solid line) and odd modes (dashed line) are shown for fixed values of 
kyHd = 7r,7r/2 and tt/4 (top, middle and bottom, respectively). Curves of constant 
effective Richardson number (dotted line) are overlayed. Riejff = const, is a good 
stability proxy at low S^/Sg. 




result of Miles (1961) and Howard (1961). At larger E^/E^, the layer is too heavy to be easily 
disturbed by a global mode. In this case, with /i S> 1, the stability limit (for either even- 
or odd-symmetry modes) is Rieif = kyH^/fi. So, at given Rie^f and Hd/Hg, the increasingly 
short wavelengths that are required for instability with increasing E^/E^ would make the 
disturbance only a surface phenomenon. Global instabilities only become possible when 
Hfi/Hg, and hence Rief^, decreases. Garaud & Lin (2004) find, similarly, that the critical 
value of Ri decreases as E^/Eg increases above ~ 0.01, for the value of vo,max/cs we adopt. 

For large E^/E^ and fixed kyH^, Figure 2 shows that the maximum unstable Hd/Hg 
has the same wavelength for both the even and odd modes. But the growth rates are not 
the same, as shown in Figure 3. As the layer thickness decreases at constant surface density, 
even modes at a given kyH^i will develop more rapidly than the corresponding odd modes. 
Or, conversely, an even mode with longer wavelength (and more efficient dilution of the 
dust layer) will become excited at the same time as a (less-efficient) odd mode with shorter 
wavelength. This can be expected, since the incomprcssibility condition implies, for the odd 
mode, that a horizontal converging flow must be generated in the midplane at the position 



Q..0,Q1 0.-010 0,tOQ.. 1..,Q.aQ 

/ 

Fig. 3. — Instability growth rates for kyHd = tt/4 and Qp = 0. Lines show the 
even (sohd) and odd (dashed) modes for vahies of the growth rate (down from the 
top) iOj/Q,K = (stability edge), 1 and 10. Dotted lines follow = const. 



where the interface displacement makes the layer thicker {5z/z > 0). Since the even mode 
docs not need to channel energy into this requirement, it can grow faster. On the other 
hand, since the odd mode does not need to vertically displace the whole layer, it can remain 
unstable for larger Ha/Hg values, at which gravity is stronger at the layer interface. 

3.5. Even Mode, Including Coriolis Forces {Vtp ^ 

In §3.2 and §3.3 we studied the solutions of equations (23) and (30) obtained when the 
Coriolis forces arc neglected. While the solutions presented were for the exact equations, we 
have found that taking the long- wavelength limit kyHd <^ 1 yields a good approximation to 
the stability even at values of kyH^ relatively large (~ 10). When the Coriolis terms are 
included {Vlp ^ 0), the dispersion relation retains its simple character when kyH^ <^ 1, so 
for this section we concentrate on that limit. ^ 



•^For short wavelengths, kyHd ^ 1 ^ iSLn\i{K2Hd) ~ 1, and equation (23) yields the same dispersion 
relation as the two layer problem with rotation studied by Chandrasekhar (1981, §105); see also Huppert 
(1968). 
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Fig. 4. — Roots of equation (32), as a function of Hd/Hg, for kyH^ = tt/3 with 
/i = 1 (a) and /i = 10 (&). Filled circles show 3?(w), and open circles show 3(c<j). 
Vertical lines mark critical Hd/Hg values (see Appendix). There is no solution for 
Rie^ > {kyHdf, i. e. large Hd/Hg. 



In the long wavelength regime, equation (23) can be written as 

/ 4Q2 \ V2 

(1 + /i) (a; - kyVfkyHd + a;M 1 - ^ j = kyg^l, (32) 

where we require the real part of the square root to be positive. The procedure used to 
solve equation (32) is presented in the Appendix, and the corresponding roots are shown in 
Figure 4. The number and nature of the roots can change only at critical Hd/Hg values (as 
described in the Appendix) , marked in Figure 4 by the vertical lines. 

For the parameters presented in Figure 4b, there is a range of values of Hd/Hg with 
only real solutions, which implies stability. Nevertheless, we consider the existence of this 
gap of limited practical significance since the gap closes for longer wavelengths (which yield 
more efficient vertical mixing). More significant is the fact that there are no solutions, real 
or complex, when (kyV/flKy < /^/(l + Z^), or kyVo^max/^K < (l + Vy")^^^ (right hand side of 
Figure 4). The transition critical point corresponds to uj/Qk = 0. For given (order unity) 
kyHd, and any I^d/^g, sufficiently small /j, guarantees that this condition holds. Thus, using 
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equation (26), a necessary (but not sufficient) condition for instability at long wavelengths 
appears to be 



Comparison of equation (33) and its equivalent for the Qp = case (cq. [28]) shows that a 
larger Rie^f is required for stability when the Coriolis terms are turned on: for kyHd = n/S, 
the critical value for Rieif is 1.1 when Qp ^ 0, compared to 0.32, 0.14, and 0.09 for fi = 1,5, 
and 10, respectively when = 0. Notice that both conditions reduce to the same criterion 
for instability when kyHd and <^ 1. However, small and Riejj < {kyHd)"^ < 1 are satisfied 
only for T^d/'^g ^ 0.01, i. e., for cases with much lower metallicity than that expected in the 
proto-solar nebula. 

An interesting feature of equation (33) is that (unlike eq. [28]) it does not explicitly 
depend on the dust abundance, ii. This is due to the fact that, at the onset of instability for 
the Qp ^ case, the phase velocity vanishes [3?(a;) — > 0]. This has two consequences. First, a 
larger fraction of the heavier layer's momentum is available to help destabihze the fiow (and 
so, the layer remains unstable for larger Kieff). And second, a vanishingly small lv means a 
small Ki (eq. [20]), so that vertical pressure gradients in the gas layer are very small. So, as 
the dust/gas mixture flows through the sinuous focus of the perturbed midplane layer, only 
the centrifugal force in the vertical direction {V'^ky 6z) and gravity {Qj^ 6z) are involved in 
the force balacing, while gas pressure gradients are neglegible. In contrast, for the Qp = 
case, 3?(a;) is non-vanishing at the onset of instability, yielding a finite Ki value. This implies 
non-neglegible vertical pressure gradients, which introduce a dependence on the dust-to-gas 
density ratio in the force balance. Small or vanishing 3?(a;) is possible only in the D,f ^ 
case, when the Coriolis forces (from vi^) are able to partly of fully compensate for azimuthal 
pressure forces. 

We can alternatively write the instability condition for Jl^^ 7^ in terms of more basic 
physical parameters by considering equation (26) in the large and small /i limits: 



Rie^ < {kyHdf. 



(33) 




(34) 



Notice that the minimum stable layer thickness approaches a constant value as J^d/^g in- 
creases, as shown by the dotted lines of Figure 2. This stands in contrast with the instability 
criterion for the flp = case (see equation [25] and Figure 2), for which {Hd/ Hg)crit oc 
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(Ed/Ep) as the surface density increases. Thus, for large the instabihty criterion is 
< vo^rnaxky for 0, Compared to {/ikyHdY/'^flK < VQ,maxky for Q,f = 0. 

While we have concentrated on the long- wavelength hmit in the discussion of this section, 
we note that the general dispersion relation of equation (23) can easily be solved in the limit 
|a;| — > 0. When 7^ 0, an unstable solution will appear when Rig^ passes below the value 



Rw = {kyHdYtscnh 



Vi-4Qy(W' 



kyHd 



(35) 



with growth rate oc (Ricrit — Rieif)^F- Thus, the necessary criterion for instability in the 
long-wavelength limit is quite similar to a more general necessary and sufficient criterion for 
a marginally- unstable mode to appear. Since tanh(x)/x > 1, BlcHt is in practice slightly 
larger than the value deduced in the long- wavelength limit. 

We note that for the odd-symmetry case, the critical Rie^f is given by the same expression 
as above, with coth(x) substituted for tanh(x); for long wavelengths, the instability criterion 
is therefore {kyVY > Aflj, + klgHdid,/{l + 11), or, using ~ 



This is always smaller than than the criterion for the even-symmetry case, and Ricrit — ^ 
{kyHd)y[4+{kyHdf] for/x> 1. 

Finally, we note that the conclusions based on the odd-symmetry three-layer system 
are similar to those for the two-layer system considered by Chandrasekhar (1981). For the 
two-layer system with rotation, instability first becomes possible when 



41^^ + ky 



l + A* 



1/2 



(37) 



This can also be written, taking Vtp fix, as 



Rie^ < /^^""^^'^ (38) 

^ 2(1 + /.) + V4(l + /.)2 + (A;,i/,)V 

For /X » 1, this yields Ricru — {kyH^y /[2+ -y/^+Jk^Hap], very similar to the odd-symmetry 
three-layer case. 
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4. NUMERICAL SIMULATIONS 

4.1. Basic Numerical Methods 

With physical understanding developed via the instability analysis of the simple three- 
layer model in §3, we now turn to fully non-linear simulations with a continuous dust distribu- 
tion. Our solutions of equations (2)- (5) are obtained using a version of ZEUS (Stone & Nor- 
man 1992), a finite difference, time explicit, operator split, eulerian code for numerical hydro- 
dynamics. For the purpose of this investigation, we implemented source terms corresponding 
to the Coriolis force and the background radial pressure gradient. We treat this background 
radial pressure gradient as constant in time and space, —dPo/dR = 2flFPgoVQ,max- 

As described in §2, in the strong-coupling limit both the gas and dust separately obey 
the continuity equation (cqs. [2] and [6]). We implement this requirement in ZEUS by 
updating using the same transport algorithm as used for the total density, p = pg + p^. 
We have tested this implementation by advecting a one dimensional, gaussian-shaped "dust 
cloud" in pressure equilibrium with its surroundings. A constant velocity was given to the 
full grid (200 zones) which advected the dust along five cloud-lengths, and across the periodic 
boundaries. This experiment was repeated for each coordinate direction of the code. In all 
cases, the total dust mass remained constant and the rms deviation from the initial gaussian 
shape was only 0.4%. 

Since ZEUS was not designed to solve incompressible fiow problems, we need to verify 
that its algorithms yield acceptable solutions in the regime of interest, and in particular, can 
properly evolve Kelvin-Helmholtz instabilities. To test this, we simulated a velocity shear 
slab with a discrete jump in velocity and density, similar to the setup in §3, with the Coriolis 
terms turned off. These simulations were performed in a two-dimensional {y — z) plane with 
initial velocity v = Vy{z)ey. The initial conditions were perturbed with an eigenfunction 
of the even mode, having a specified wavelength. We then measured the growth rate of 
the specified perturbed mode, while it remained in the linear phase and before other modes 
(faster growing with smaller wavelengths, seeded by grid-size noise) , interfered in any obvious 
way. We compared this growth rate with the analytical value for a variety of perturbation 
amplitudes, obtaining an agreement of the order of 10%. 

We now move on to the initialization of our full non-linear, continuous simulations. 
For each model, we set up a disk atmosphere in vertical hydrostatic equilibrium, and with 
azimuthal velocity Vy{z) compensating radial gravity and the background pressure gradient. 
As described in §2, setting the value of ii{z) — pd{z) / Pg{z) defines all the hydrodynamic 
variables in the equilibrium state. 
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4.2. Experiments With Prescribed fj,{z) 

The experiments described in this section are initiahzed by simply defining n as an 
arbitrary function of z, namely, 

IJ,{z) = iio sech^ {z/ Hd) , (39) 

with Ho = /x(0). Three parameters uniquely define a model: /iq, vo,max/cs and Hd/Hg (for 
consistency with §3, in practice we shall use ^df^g instead of Hq). In these experiments, 
the simulation domain is a square grid of size 20Hd (257 x 257 grid points), with periodic x- 
and y-boundaries, and closed 2;-boundaries.^ The equilibrium state is perturbed by random 
velocities with a maximum amplitude of 10~^Cs. In all the cases, vo^max — O-lc^ which, for 
MSN parameters, corresponds to a radial location at R 3.5 AU. 

Figure 5 shows results from the simulations with Hd — O.OlHg and = 0.008Eg, both 
with the Coriolis terms turned off (left column; Qp — 0) and with the Coriolis terms turned 

on (right column; Qp ^ 0). The density snapshots presented arc t = 1.5,2.2 and 2.9torb 
[Qp = 0), and t = 1.8, 2.3 and 2.8torb (^f 7^ 0), where torb = ^tt/Qk is the orbital time. The 
initial linear growth phase is very similar for both cases, with a sinusoidal displacement of the 
dust layer consistent with the even (midplanc-crossing) mode of the instability. Notice that 
the instability develops three wavelengths in the computation domain, which corresponds to 
kyHd ~ 0.95. We also performed simulations with 3 and 5 times larger grids in the y direction 
(at lower resolution), and found that the same approximate kyHd mode predominated. The 
fastest growing mode in the analysis presented in Ishitsu & Sekiya (2002, 2003) has a similar 
kyHd value. 

In the D,F — case, once the displacement is large enough, drag with the dust-poor 
fluid erodes the tips of the displaced layer, generating famihar backward-facing KH rolls. 
The flow around the rolls pushes the dust toward the peaks, where it accumulates until 
the layer breaks. After the wave breaks, the dust mixes with the surrounding gas, lowering 
the dust abundance and the velocity shear, and generating a more-or-less homogeneous layer 
with a thickness similar to the wavelength of the linear growth. (For some parameter choices, 
this layer is in turn unstable, and again developes a linear growth phase.) Remnant non- 
organized motions generate difi'usion that further increase the layer thickness over longer 
time-scales. 

In contrast to the ftp = ^ model, after the linear growth phase the flp ^ case has a 



"^We also tested open z-boundaries, and found no significant differences. 
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Fig. 5. — Sample snapshots from a simulation with prescribed fJ.{z), both with 
Coriolis terms turned off (Oj? = 0, left column) and turned on (ilp / 0, right 
column). The grayscale shows dust density and the arrows show the velocity in the 
dust-free reference frame; the arrow in the box represents a magnitude of fo-maa / Cs ■ 
Coordinate axes are in units Hg. Both models have = O.OlHg, 11^1 = O.OOSSg, 
and VQ^rnax = Clc^. Snapshots correspond to t = 1.5,2.2 and 2.9 tcrb for the left 
column, and t = 1.8,2.3 and 2.8 tgrb for the right column. 
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Fig. 6. — Snapshots from the 0,f ^ simulation shown in Figure 5, at t = 1.8 
(left) and 2Morb (right). The grayscale shows dust density, and lines show contours 
of Vx'- thin solid for > 0, thick solid for Vx = 0, and dashed for Vx < 0. As the 
perturbation grows, the dust layer experiences the slower-moving high latitude gas 
as a head-wind. This gas steals angular momentum from the dust layer, which flows 
radially inward, while the dust-poor gas flows outward. 



period of very rapid evolution. Unlike the = Q case, the dust accumulates in the midplane 
at the nodes of the displaced layer. As the vertical displacement grows, the leading faces 
of the perturbed dust layer encounters the slower-moving high-latitude gas. As the gas is 
pushed by the dust, it gains angular momentum and flows radially outward (see Figure 
6). Correspondingly, the dusty layer looses angular momentum and flows radially inward, 
increasing its azimuthal velocity. As a consequence, the displaced peaks move faster than 
the gas in the midplane, generating forward leaning structures. When the initial density 
layer looses integrity, it forms a new layer thicker than the one generated in the 0,^ — case. 
The subsequent diffusive growth also develops more rapidly. 

Nearly every set of parameters we tested triggered growth dominated by an even mode, 
consistent with expectations from the discussion in §3.4. In order to induce the development 
of the odd mode, we forced the vertical symmetry by placing a boundary of the simulation 
at 2; = 0, for the Qp = case, at the same spatial resolution. Results for this experiment 
are shown in Figure 7 (the grid is doubled in this plot for easier comparison with Figure 5). 
There are several differences evident between these examples of even and odd modes. First, 
the odd mode dcvelopcs longer wavelengths, fitting only two wavelengths in the domain, 
evolving later to a single wavelength. Second, despite the longer wavelength, the extent of 
the layer after the non-linear growth phase is thinner, with the dust accumulated at the 
midplane. Third, this (and other) odd mode simulations evolve much more slowly than the 



Fig. 7. — Simulation with prescribed fi{z) for the i}p = case, forcing the odd mode 
of the instabihty by placing a reflecting boundary at z = 0. From top to bottom, 
snapshots show dust density (grayscale) and velocity field (vectors) at t = 2.5, 5.0 
and 7.5 torb- Notice the smaller vertical extention of the grid, compared with Figure 
5. 



even mode ones. This odd mode simulation has a much larger latency time, despite the fact 
that the initial perturbations are a factor of 10 larger than the even mode example. The 
dominant wavelengths that develop for the even modes have larger growth rates than the 
ones for the odd modes: for this case, 0.5SQ,k vs. 0.05^/^ (the procedure to measure these 
growth rates is described in §4.3). 

In addition to the sample cases shown, we have performed a large number of similar 
models with both flp = and flp 0, and either full or "half" grids. We have covered a 
range H^/Hg = 0.001 through 0.08, with a constant surface density ratio '^d/^g — 0.008. 
Overall, we find a confirmation of the above results: the layer evolves more rapidly when 
fip 7^ 0, and more slowly when we force the development of the odd mode. We also found 
that, for Vtp = and this E^/E^, value, the layer becomes stable for Hd/Hg > 0.02. When 
flp ^ 0, the layer appeared to be always unstable, although the evolution times are quite 
small for the thickest layers. 
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4.3. Experiments With Ri = const 

The initialization procedure adopted in §4.2 is quite simple, but also quite arbitrary in 
terms of the detailed functional form of the dust distribution. As an alternative approach, 
we can adopt a prescription for the initial conditions which satisfies additional physically- 
motivated constraints. A natural choice is to set Ri to a constant value everywhere in the 
flow; this choice is motivated by the common conception that low-Ri flows are unstable, and 
high-Ri flows are stable. In particular, Sekiya (1998) and Youdin & Shu (2002) have argued 
that the disk may evolve to a state where Ri = 1/4 (see §1). 

Using the definition of the Richardson number in equation (1), and equations (10) and 
(15), we can obtain a differential equation for /i (for a given constant Ri): 



d// _ -z{l + nf 

dz 2Ri H^{vo^rnax/Csy 



1 + Wi + 



4Ri(^;o,^a^/c^)^ 
l + A* 



(40) 



An approximate solution to this equation can be found provided that 4Ri (i»o,r„aa;/cs)^/(l + 
/x) -C 1 (true for the parameter space of interest). Using this approximation. 



l + HoJ ^i{Vo,max/Csy 



-1/2 



(41) 



Again, three parameters define a model: /iq (or alternatively, S(i/Sg), VQ^max/cg, and Ri. This 
equation sets a finite extension for the dust layer. Setting ij,{z) = 0, we obtain 
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j^Ri 


^0,max 











1 - 



;i+/^o)^ 



1/2 



(42) 



We use this setup to study the stability of the dust layer under varying conditions. 
Cases with and without the Coriolis terms were performed. For our array of simulations, 
we adopt VQ^max/cs — 0.1, and explore a range of initial values of S^/E^ and Ri. Our 
parameter grid covered E^/Eg = 0.006 to 0.178, and Ri = 0.1 to 1.0 for the Q,f — ^ cases, 
Ri = 1.0 to 10.0 for the VLp ^ cases. The domain of these simulations extended ±lQzmax 
in and -^2zmax in z (we did not need a large z extension since these experiments were 
designed to follow only the linear part of the instability, as opposed to the full evolution of 
§4.2). Our standard resolution was using 200 x 40 grid points, with higher resolution models 
tested for selected cases in order to verify the results. Each simulation was followed through 
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Fig. 8. — Left: Fourier transform of Vz along the y direction, as a function of z, 
for the Ri = 0.158, Ti^/Tig = 0.056, fip = case, after 4.27 orbital times. Right: 
Amplitude of the A = 0.26-ffg mode, in arbitrary units, as a function of time, for 
2; = 0. The linear growth phase of the instability (solid line) shows a growth rate 

uji = o.unK. 



t/torb = 4:00zmax/ Hg, whlch coriesponds to values between 9 and AOtorh for Ri between 0.1 
and 1.0. The initial equilibrium setup was perturbed with random velocities with a maximum 
amplitude of lO'^c^. 

After exploring several possibilities, we found that a practical way of assessing the 
stability of the dust layer was by calculating the Fourier transform of Vz along the y direction, 
as a function of z. When the layer is unstable, the fastest growing mode is easily identified 
from this plot; its amplitude can then be examined as a function of time. In most cases, the 
linear growth phase of the instability is clearly evident, and a growth rate can be measured. 
An example is presented in Figure 8. 

As expected, the empirically-determined region of stability in the parameter space does 
not have a sharp edge, but instead growth rates decrease until our procedure becomes too 
insensitive to measure them. For practical purposes, we denote as "stable" those simulations 
that do not show growth within the allotted run time, or else for which the measured growth 
rate is lower than a fiducial value, set to uj — O.OSJlx; this growth rate corresponds to a 
ten-fold increase in amplitude over 10 orbital times. 

In Figure 9 we present results for an array of values in the E^/Eg — Ri parameter space. 
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Fig. 9. — Stability in the dust layer for a matrix of values of the dust surface 
density Sd/S^ and Richardson number values Ri, (a) with Coriolis terms turned off 
{p,F = G) and (6) with Coriolis terms turned on {ytp ^ Q). Stable cases, based on 
simulation growth rates, are marked with an "S" and unstable cases with a "U." 
Models with a measured growth rate toj < 0.05^2 ft- are considered stable. The lines 
across these plots trace the value of the midplanc dust fraction, (from the left) 
ji = 0.3, 1.0 and 10. Note the different ranges of Ri in (a) and (6). 



noting whether each model is stable (S) or unstable (U) according to the above criteria. 
Generally speaking, the Vtp = ^ case shows stability for Ri > 0.3 (Figure 9a). Similarly to 
§4.2, we also forced the development of the odd mode of the instability by placing a reflecting 
boundary at 2; = 0. For the odd modes, the edge of stability also appears consistent with the 
criterion Ri = 1/4, although the measured growth rates are much lower than for the even 
mode. 

The main difference between the Vtp = {) and Vtp ^ {) cases (Figure 9b) is that stability 
in the Vtp ^ {) case requires a much larger value of the Richardson number, near Ri ~ 5. It 
is important to stress that this threshold value for Ri depends on the value chosen for the 
fiducial minimum "unstable" growth rate, since there is some subjectivity in the procedure 
to measure growth rates and the slope of uji as function of Ri is quite shallow when 7^ 
(see Figure 10). Nevertheless, the layer is quite clearly unstable for values of Ri significatively 
larger that 1/4. Since Ri can be thought as a proxy for the layer thickness (eq. [42]), the 
condition that a larger Ri is required for stability when Q^p ^Qis consistent with the results 
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Fig. 10. — Comparison of growth rates for the Up = (dashed hnes) and Up 7^ 
(soUd hnes) cases, in units of fix- The growth rates are generally higher when 
Coriolis terms are turned on. 



of §4.2, in which we found that the final dust distribution is thicker when the Coriolis terms 
are turned on. 

As a coda for this section, we comment on the constant Ri distribution adopted here. 
Sekiya (1998) studied the profile determined by a constant-Ri dust distribution and noticed 
that the space density develops an infinite cusp if the S^/Sp value is large enough. Youdin 
& Shu (2002) interpreted this result as a hint that the gas distribution can only sustain 

the weight of a mass of dust similar to its own mass, and if surpassed, the dust would 
precipitate. By numerically integrating our equation (41) for a given Ri = const., together 
with the hydrostatic equilibrium condition in equation (15), we find that 



for yUo ^ 1 (see also eq. [22] in Sekiya 1998). This means that, for Ricrit of order unity, 
Tin/Tig > Vo^max/cs implies an exponential growth of the midplane dust abundance. The 
onset of a cusp in density would therefore occur (under the condition Ri = const.), for a 
range of radii ~ 1 — lOAU (see eq. [14]), when T,d/T,g exceeds ~ 0.1, i. e., for roughly an 




(43) 
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order of magnitude dust enhancement. However, this conclusion appears to be an artifact 
of the Ri = const, condition. If Ri is permitted to vary with height, with Ri becoming large 
near the midplane, then /xq varies only linearly rather than exponentially with T^dfEg, and 
a cusp need not form. 

5. SUMMARY AND CONCLUSIONS 

In this work, we have studied the stability of gaseous disks with vertically varying dust 
abundances and accompanying vertical velocity gradients. We present both linear analysis 
of a simple discrete three-layer distribution, and full numerical simulations of a pair of 
continuous distributions. Wc focus on identifying and intercomparing criteria for instability 
when the Coriolis force terms arc considered (Qp ^ 0) or disregarded {Qp = 0). Our 
results indicate that the Richardson number Ri is still a good general discriminant for the 
stability /instability boundary when Qp ^ 0, but the critical value is larger than the classical 
result Ri = 1/4 for Qp = 0. 

As with the well known two-layer problem, the critical Ri for the discrete thrcc-Iaycr 
distribution depends on the wavelength of interest. We argue that wavelengths of order of 
a few times the layer thickness {kyH^ ^ 1) are the most relevant, since smaller wavelengths 
will yield less dust mixing in the non-linear regime, and larger ones have smaller growth 
rates. In addition, for realistic vertical distributions with continuous variations of Vy{z) and 
^i{z), waves with kyHd ^ 1 would not be unstable. For kyHd = n/S, RUff < 1-1 is a 
necessary (but not sufficient) condition for instability for the Qp ^ case. In general, when 
Qp 0, marginally-unstable modes first appear when Rie^f drops below {kyHd)'^. Since Rigjff 
is a proxy for the layer thickness, this means that, when Coriolis forces are considered, the 
dust layer is thicker than in the Vtp = case at the onset of instability. At S^/Ep = 0.01, 
the increase in thickness is ~ 50%, whereas for Sd/^^s = 0.1, the increase is nearly a factor 
of three. For the VQ^rnax/cs — 0.1 value adopted throughout this work, R — 3.5 AU (for 
MSN parameters), and GI will set in only when ji > 187 (Sekiya 1998). The thicker layer 
at the "edge" of KHI imphes that the dust-to-gas surface density ratio necessary for GI 
(neglecting vertical self-gravity) would increase from the already high value of S^/E^ = 1.1 
(eq. [25]) to S^/Sg = 15 (eq. [34]), when Coriolis effects are included.^ While we certainly 
do not believe that our results (assuming full coupling) can be accurately extrapolated into 



^For the continuous dust distribution studied by Garaud & Lin (2004), the surface density needed for GI, 
with these parameters and neglecting Coriohs forces, would be S^/E^ « 0.4, only slightly lower than the 
value 1.1 for our 3 layer model with Qp = 0. 
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the dust-dominated regime, this comparison gives a sense of the strong change that can be 
expected. The astronomical imphcation is, of course, that GI is even more difficult to achieve 
than previously thought. 

We note that consideration of even-symmetry (in Vz, i. e. midplane-crossing) modes is 
important for this problem, although some previous analytic studies (Ishitsu & Sekiya 2002, 
2003, for example) have focused only on odd-symmetry modes. The even-symmetry case 
becomes unstable at a larger value of Rigjff than the odd-symmetry case [which is similar to 
the two-layer system analyzed by Chandrasekhar (1981)]. 

The results from our simulations are consistent with the main conclusions drawn from 
the simple three-layer stability analysis. Based on an array of simulations performed that 
were initialized with Ri = const, vertically, we found that Ri ~ 5 is the critical value found 
for the Qp ^ ^ case. This critical Ri value should be considered an approximate result, 
since the procedure used to evaluate stability from the simulations involves a certain level of 
subjectivity. Nevertheless, the dust layer is clearly unstable for Ri values significantly larger 
than 1/4. As with the discrete case, the larger Ri value also imphes a thicker dust layer; for 
Erf/Sp = 0.01 at the threshold of instability for Qp ^ 0, the height containing half the dust 
surface density increases by a factor ~ 3 compared to the corresponding non-rotating model. 
Of course, both this and other results based on the strong-coupling assumption should be 
considered carefully when ^ 1, i. e. for large E^/Sp. 

Youdin & Shu (2002) point out that precipitation of solids is hkely to occur in layers 
for which /i » 1. For vertical profiles with Ri = const., this situation arises first at the 
midplane, where a density cusp forms for E(i/Sg exceeding some critical value (of order 
Vo^max/cs)- They therefore conclude that, for Gl to develop, E(i/Sg needs to be increased 
from the MSN value by a smaller amount than estimates without a density cusp. Since our 
results imply that stability fails at larger Ri (and hence a thicker dust layer) than the value 
assumed by Youdin & Shu (2002), the total dust enhacement S^/Sg required before a cusp 
is predicted (and precipitation develops) would nominally be larger than they concluded. 
However, it is not clear that the condition for cusp formation can in fact be obtained and 
written in a simple way in terms of E^i/Eg, for the reason we now discuss. 

States which have Ri = const, have been considered of particular interest because they 
are everywhere at a nominal threshold for instability. However, the proposal that profiles 
should evolve to such states implicitly depends on an assumption that the localized failure of 
the KH stabihty criterion would lead to locahzed turbulence. While we chose a constant-Ri 
setup in the models of §4.3 to facilitate comparison with previous work, the results of the 
simulations in fact call this implicit assumption into question. When exploring Ri(2;) for 
the simulations in §4.2, in some of our models the Ri > 1/4 stability condition failed only 
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for \z\ < 2Hd. Nevertheless, the fastest growing mode in these models had A ^ 7Hd, and 
the non-linear evolution involved the whole layer, yielding dust mixing at heights over 10 
times the initial value. Furthermore, when the simulations were followed until the dust 
distribution was somewhat homogeneous in the azimuthal direction, the vertical profile did 
not appear to have a constant Ri value. 

It is our view that, as exemplified by these models, a local failure to satisfy a stability 
criterion which is stated in terms of local quantities need not trigger only local instabilities, 
i. e. of wavelengths similar to the thickness of the "failing" region. Instead, the wavelength 
will be similar to the overall velocity gradient, and the height of the waves when they 
break will be similar to the wavelength. The re-adjustment of the dust distribution that 
is generated will thus be global, including regions that initially have Ri > Ricru- It is, 
therefore, not obvious that the "final," or quasi-equilibrium, distribution should be accurately 
described by a local quantity. It may be that the marginal state that develops as a result 
of the competition between dust settling and shear-driven turbulence in fact has a relatively 
uniform distribution of dust in a central layer, rather than a cusp out of which solids can 
precipitate. Of course, if particles grow large enough prior to the onset of global dust-layer 
turbulence, then they may continue to settle faster than the growth rate of the instabilities. 
Since particle growth depends on uncertain sticking efficiency, however, this question remains 
open. Observations with newly deployed infrared telescopes should shed some light on the 
true vertical distribution of solids in proto-planetary disks. 

For the technical modeling in this work, we have adopted a number of simplifications 
and idealizations. While we believe that our main conclusions are not sensitive to these 
assumptions, it is appropiate to review what they are. One such simplification is that we 
consider only models with — 0, and neglect radial shear and tidal forces. We focused on 
kx — because these modes typically have the fastest growth rates and largest instabihty 
thresholds in terms of Ri. Had we included shear, then k^ would have grown in time as 
dkx/dt = (3/2)f2xfcy, which would ultimately limit growth of perturbations of a given ky. 
In future work, it will be important to examine how, for example, the total amplification of 
shearing KH wavelets depends on Ri. 

Other idealizations include neglecting the slip of gas relative to dust, treating the large- 
scale radial gradient of pressure as fixed in time, and neglecting other sources of turbulence. 
Gas-dust slip can itself lead to streaming instabilities (Youdin & Goodman 2005). However, 
either with or without full coupling of the gas and dust, the instabilities that develop depend 
on the gas having consistent sub-keplerian azimuthal velocities. This relies on having a 
consistent radial pressure gradient; in principle, this condition could be invalidated if there 
are sufficient pressure perturbations induced by turbulence or large-scale waves in the disk. 
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The ratio of such terms, ~ {SP R)/{PoX), could be significant even for moderate amplitude 
perturbations with A ~ Hg. In future work, we intend to explore these and other dynamical 
processes that may affect the settling and mixing of solids in dusty disks. 
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A. SOLUTION TO THE DISPERSION RELATION, FOR ^ 

The dispersion relation for the even-symmetry mode, equation (23), in the long- wave- 
length regime, is given in equation (32). To investigate the nature of the solutions of this 
non-algebraic equation, we follow a graphical procedure similar to that in Chandrasekhar 
(1981, p. 498). Define 
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In these variables, the dispersion relation reads 



l + li 



By simultaneously considering the equation 
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Fig. 11. — Locus of the real solutions to equations (A4) (thick solid line) and (A5) 
(thin solid line), for kyH^ = 7r/3 and /x = 10. Straight lines C ~ ^7 = const, are 
related to critical points of the problem (see text). For given values of kyH^, /x and 
/cs, varying values of ^ - 77 correspond to different values for Ha/Hg. 



together with the relation |^ — r/l = kyV/flK = const., we can obtain a fourth order poly- 
nomial instead of a non-algebraic equation. All the solutions of the polynomial must be 
solutions either of equation (A4) or (A5). 

The thick solid line in Figure 11 shows the locus of the real solutions of equation (A4), 
which we shall call the P-branch, while the thin solid line shows the real solutions of equation 
(A5), which we shall call the ^'-branch. For clarity, we show only the ^ > side of the 
diagram. In addition. Figure 11 shows several lines representing = const, which relate 

to a series of critical points (^, 77) for the problem as follows: 

1. The ends of the P-branch he at (±/3, ±770); these are the intersections with the dash- 
dotted lines ^ — ri — P ±rio in Figure 11. Here, tIq — /x/(l -|- //). 



2. The maximum extention of |.^| for the P-branch lies at (±,^05 0); this is the intersection 
with the long-dashed line ^—r] = ^o in Figure 11. Here, = \/ {nkyHdY + /?^/4-F/3^/2. 
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3. The point of the P-branch with the maximum value of |^ — 77I, {Ct,Vt)', this is at the 
intersection with the dotted hne in Figure 11. 

4. The (^,77) = (0, — ?7o) point, where r]Q = + /x); the short-dashed hne i — i] = —rjQ 
in Figure 11 runs through this point. 

Any of the four solutions to the equivalent fourth-order polynomial must be a solution 
of either equation (A4) or (A5), and the number or nature of the roots of each equation can 
change only at the critical points outlined above. The roots corresponding to equation (A4) 
are shown in Figure 4, for an array of Hfi/ Hg (or equivalently |^ — ?7|) values, given /i — 1 
(a), and /i — 10 (b), and kyHa — tt/S. 

On the left the diagrams in Figure 4, at small H^/Hg (or on the right of Figure 11, at 
large |^ — there is a pair of complex roots. In Figure 4b, these two complex roots become 
real as — r]\ becomes less than — Vt since then a — r^l = const, line crosses the P- and 
S'-branches four times, implying all four roots arc real. Moving further right in Figure 4b 
(or left in Figure 11), one of those roots becomes complex when — ^^l < (3 + rjo, and the 
other real root disappears when \i — ri\ < P — rjo. In Figure 4a, at lower n, a complex root 
is present whenever a real root exists. At the far right in Figures 4a and 4b (corresponding 
to the region of Figure 11 with |^ — r^l < r^o), there are no solutions. 

There is a range of values of |^ — ?7| with all real solutions (implying stability for that 
set of parameters) if the P-branch has a section with |dr]/d{| < 1. When that happens (as 
for the fjL — 10 case of Figure 4b), the points of the P-branch with |d7;/d^| = 1 he very close 
to the (/3, — ?7o) and (^0, 0) points. So, a good approximation for the opening of such stability 
gap is ^0 > P + Vo- While it is straightforward to find a precise criterion for the existence 
of this gap, it is of limited practical significance since the gap closes for longer wavelengths 
(which yield more efficient vertical mixing). We consider more significant the fact that there 
are no solutions, real or complex, for |^ — < fjo, since this yields a necessary (but not 
sufficient) condition for instability. 
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